#include "Box3.hpp"
#include <limits>
#include <map>
#include "MeshAlgorithm.hpp"
#include <queue>
#include <iostream>
#include "EdgeMesh.hpp"
#include "PolyMesh.hpp"
#include "CreaseStructure.hpp"
#include "DebugHelper.hpp"
#include "MeshAdjacency.hpp"
using namespace std;
using namespace MeshAlgorithm;
using namespace Algorithm;
//
//void MeshAlgorithm::BuildEdgeGraph(EdgeGraph&              aGraph,
//                                   const std::vector<Edge> aEdges,
//                                   int                     aVertexCount)
//{
//   for (size_t i = 0; i < aEdges.size(); ++i)
//   {
//      add_edge(aEdges[i].first, aEdges[i].second, aGraph);
//   }
//}
//struct CycleVisitor : public boost::default_dfs_visitor
//{
//   //typedef MeshAlgorithm::EdgeGraph Graph;
//   //typedef Graph::edge_descriptor Edge;
//   //typedef Graph::vertex_descriptor Vertex;
//   //void initialize_vertex(Vertex&, Graph& g){}
//   //void start_vertex(Vertex&, Graph& g){}
//   //void discover_vertex(Vertex&, Graph& g){}
//   //void examine_edge(Edge& e, Graph& g){}
//   //void tree_edge(Edge& e, Graph& g){}
//   //void back_edge(Edge& e, Graph& g){}	
//   //void forward_or_cross_edge(Edge& e, Graph& g){}
//   //void finish_vertex(Vertex&, Graph& g){}
//};

struct CycleGraph
{
   typedef std::multimap<int, int> EdgeMap;
   std::set<int> unvisited;
   std::multimap<int, int> adj;
};
void ClearLinked(CycleGraph& g,
                 int         aV)
{
   std::vector<int> verts;
   verts.push_back(aV);
   int prev=-1;
   while (! verts.empty())
   {
      int v = verts.back();
      g.unvisited.erase(v);
      verts.pop_back();
      CycleGraph::EdgeMap::iterator i = g.adj.lower_bound(v);
      CycleGraph::EdgeMap::iterator endI = g.adj.upper_bound(v);
      for (; i != endI; ++i)
      {
         if (i->second == prev) continue; // quick culling prev
         if (g.unvisited.find(i->second) != g.unvisited.end())
         {
            verts.push_back(i->second);
         }
      }
      prev = v;
   }
}
bool CycleDfs(CycleGraph&       g,
              int               v,
              std::vector<int>& cycle,
              int               aMaxCycle)
{
   int nextV = -1;
   bool ok = true;
   int prev(-1);
   while (ok)
   {
      std::set<int>::iterator iter = g.unvisited.find(v);
      if (iter == g.unvisited.end())
      {
         return true;
      }
      g.unvisited.erase(iter);
      CycleGraph::EdgeMap::iterator i = g.adj.lower_bound(v);
      CycleGraph::EdgeMap::iterator endI = g.adj.upper_bound(v);
      std::vector<int> branches;
      int branchCount(0);
      int branch;
      for (; i != endI; ++i)
      {
         if (i->second == prev) continue;
         if (prev == -1) { prev = i->second; continue; }
         ++branchCount;
         branch = i->second;
      }
      cycle.push_back(v);
      prev = v;
      if (branchCount == 1)
      {
         v = branch;
      }
      else
      {
         ok = false;
      }
      if ((int)cycle.size() > aMaxCycle)
      {
         ok = false;
      }
   }
   ClearLinked(g, v);
   return false;
}
void MeshAlgorithm::FindCycles(const std::vector<Edge> aEdges,
                               int                     aMaxEdgeCount,
                               std::list<EdgeCycle>&   aCycles)
{
   CycleGraph g;
   
   for (size_t i = 0; i < aEdges.size(); ++i)
   {
      Edge e = aEdges[i];
      g.adj.insert(e);
      g.adj.insert(swap_edge(e));
      g.unvisited.insert(e.first); 
      g.unvisited.insert(e.second);
   }
   while (! g.unvisited.empty())
   {
      EdgeCycle cycle;
      if (CycleDfs(g, *g.unvisited.begin(), cycle, aMaxEdgeCount))
      {
         aCycles.push_back(EdgeCycle());
         aCycles.back().swap(cycle);
      }
   }
   
}

MeshAlgorithm::VertexFaceLinks::VertexFaceLinks(int         aMaxFaces,
                                                const Mesh& m)
{
   mMaxFaces = aMaxFaces;
   mData.resize(m.verts.size() * mMaxFaces, -1);

   for (size_t i = 0; i < m.faces.size(); ++i)
   {
      const Face& f = m.faces[i];
      InsertLink(i, f.vertex[0]);
      InsertLink(i, f.vertex[1]);
      InsertLink(i, f.vertex[2]);
   }
}

struct IsDeadVert
{
   Mesh*                mMesh;
   std::vector<int>*    mUses;
   bool operator()(Vec3& v)
   {
      int offset = &v - &mMesh->verts.front();
      return (*mUses)[offset] == -1;
   }
};
void MeshAlgorithm::RemoveFacelessVertices(Mesh& aMesh)
{
   std::vector<int> useCount(aMesh.verts.size(), 0);
   for (size_t i = 0; i < aMesh.faces.size(); ++i)
   {
      Face& f = aMesh.faces[i];
      ++useCount[f.vertex[0]];
      ++useCount[f.vertex[1]];
      ++useCount[f.vertex[2]];
   }

   int vertNumber = 0;
   for (size_t i = 0; i < aMesh.verts.size(); ++i)
   {
      if (useCount[i] == 0)
      {
         useCount[i] = -1;
      }
      else
      {
         useCount[i] = vertNumber++;
      }
   }

   if (vertNumber < (int)aMesh.verts.size())
   {
      for (size_t i = 0; i < aMesh.faces.size(); ++i)
      {
         Face& f = aMesh.faces[i];
         f.vertex[0] = useCount[f.vertex[0]];
         f.vertex[1] = useCount[f.vertex[1]];
         f.vertex[2] = useCount[f.vertex[2]];
      }

      IsDeadVert pred;
      pred.mMesh = &aMesh;
      pred.mUses = &useCount;
      VertexList::iterator newEnd = std::remove_if(aMesh.verts.begin(), aMesh.verts.end(), pred);
      aMesh.verts->erase(newEnd, aMesh.verts.end());
   }
}

void MeshAlgorithm::PrintMeshInfo(Mesh& aMesh)
{
   cout << "Vertex Count: " << aMesh.verts.size() << '\n';
   cout << "Face Count: " << aMesh.faces.size() << '\n';
}

void MeshAlgorithm::SplitMaterials(Mesh&            aMesh,
                                   std::list<Mesh>& aMeshParts)
{
   std::set<int> materials;
   for (size_t i = 0; i < aMesh.faces.size(); ++i)
   {
      Face& f = aMesh.faces[i];
      materials.insert(f.materials[0]);
      materials.insert(f.materials[1]);
   }

   for (set<int>::iterator i = materials.begin(); i != materials.end(); ++i)
   {
      int mat = *i;
      aMeshParts.push_back(Mesh());
      Mesh& m = aMeshParts.back();
      for (size_t j = 0; j < aMesh.faces.size(); ++j)
      {
         Face& f = aMesh.faces[j];
         if (f.materials[0] == mat)
         {
            m.faces->push_back(f);
         }
         else if (f.materials[1] == mat)
         {
            Face flip(f);
            swap(flip.materials[0], flip.materials[1]);
            swap(flip.vertex[0], flip.vertex[2]);
            m.faces->push_back(flip);
         }
      }
      *m.verts = *aMesh.verts;
      RemoveFacelessVertices(m);
   }
}

void MeshAlgorithm::MakeManifold(Mesh& aMesh)
{
   EdgeMesh em(aMesh);
   const EdgeMesh::FaceList& faces = em.GetFace();
   const EdgeMesh::EdgeList& edges = em.GetEdges();
   std::multimap<int, int> edgeToFace;
   typedef multimap<int,int>::iterator EFIter;
   std::vector<int> edgeCounts(edges.size(), 0);
   for (size_t i = 0; i < faces.size(); ++i)
   {
      const EdgeMeshData::EdgeFace& f = faces[i];
      ++edgeCounts[f.edges[0]];
      ++edgeCounts[f.edges[1]];
      ++edgeCounts[f.edges[2]];
      edgeToFace.insert(pair<int,int>(f.edges[0], i));
      edgeToFace.insert(pair<int,int>(f.edges[1], i));
      edgeToFace.insert(pair<int,int>(f.edges[2], i));
   }
   std::vector<int> faceMarkers(faces.size(), -1);
   std::vector<int> markerCounts;
   bool doneMarking(false);
   for (int mark = 0; !doneMarking; ++mark)
   {
      markerCounts.push_back(0);
      int face = -1;
      for (size_t i = 0; i < faceMarkers.size(); ++i)
      {
         if (faceMarkers[i] == -1)
         {
            face = i;
            break;
         }
      }
      if (face != -1)
      {
         std::vector<int> neighbors;
         neighbors.push_back(face);
         do 
         {
            int current = neighbors.back();
            neighbors.pop_back();
            if (faceMarkers[current] != -1) continue;
            faceMarkers[current] = mark;
            ++markerCounts[mark];
            
            for (int ei = 0; ei < 3; ++ei)
            {
               int e = faces[current].edges[ei];
               if (edgeCounts[e] == 2)
               {
                  EFIter l = edgeToFace.lower_bound(e);
                  if (l->second == current)
                  {
                     ++l;
                  }
                  int nb = l->second;
                  if (faceMarkers[nb] == -1)
                  {
                     neighbors.push_back(nb);
                  }
               }
            }
         } while (! neighbors.empty());
      }
      else
      {
         doneMarking = true;
      }
   }

   std::vector<int> badFaces;
   for (size_t i = 0; i < edgeCounts.size(); ++i)
   {
      if (edgeCounts[i] > 2)
      {
         typedef std::pair<int, int> FaceChoice;
         std::vector<FaceChoice> choices;
         EFIter begIter = edgeToFace.lower_bound(i);
         EFIter endIter = edgeToFace.upper_bound(i);
         for (EFIter j = begIter; j != endIter; ++j)
         {
            int faceIndex = j->second;
            int marker = faceMarkers[faceIndex];
            choices.push_back(FaceChoice(markerCounts[marker], faceIndex));
         }
         std::sort(choices.begin(), choices.end(), std::greater<FaceChoice>());
         for (size_t c = 2; c < choices.size(); ++c)
         {
            badFaces.push_back(choices[c].second);
         }
      }
   }
   std::sort(badFaces.begin(), badFaces.end());
   badFaces.erase(std::unique(badFaces.begin(), badFaces.end()), badFaces.end());
   cout << "Make manifold found " << badFaces.size() << " 'bad' faces.\n";
   for (int i = badFaces.size() - 1; i >= 0; --i)
   {
      aMesh.faces->erase(aMesh.faces.begin() + badFaces[i]);
   }
   if (badFaces.size())
   {
      RemoveFacelessVertices(aMesh);
   }
}

namespace
{
   struct EdgeData
   {
      int faceCount;
   };

   //struct EdgeHash
   //{
   //   static const int bucket_size = 8;
   //   static const int min_buckets = 1000;
   //   size_t operator()(const Edge& aEdge) const
   //   {
   //      return aEdge.first | (aEdge.second << 16);
   //   }
   //   bool operator()(const Edge& a, const Edge& b) const
   //   {
   //      return a < b;
   //   }
   //};

   typedef std::map<Edge, EdgeData> EdgeMap;
   

   void InsertEdge(int a, int b, EdgeMap& aMap)
   {
      Edge e = OrderEdge(a,b);
      EdgeMap::iterator i = aMap.find(e);
      if (i == aMap.end())
      {
         EdgeData ed;
         ed.faceCount = 1;
         aMap.insert(EdgeMap::value_type(e, ed));
      }
      else
      {
         ++i->second.faceCount;
      }
   }

   typedef std::pair<Edge, int> MaterialEdge;
   typedef std::map<MaterialEdge, EdgeData> MaterialEdgeMap;

   void InsertEdge(int a, int b, int aMat, MaterialEdgeMap& aMap)
   {
      Edge e = OrderEdge(a,b);
      MaterialEdge me(e, aMat);
      MaterialEdgeMap::iterator i = aMap.find(me);
      if (i == aMap.end())
      {
         EdgeData ed;
         ed.faceCount = 1;
         aMap.insert(MaterialEdgeMap::value_type(me, ed));
      }
      else
      {
         ++i->second.faceCount;
      }
   }
}

void MeshAlgorithm::PrintEdgeInfo(Mesh& aMesh)
{
   //std::map<Edge, int> edgeIndices;
   //std::vector<Edge> edges;
   EdgeMap edges;
   for (size_t i = 0; i < aMesh.faces.size(); ++i)
   {
      Face& f = aMesh.faces[i];
      InsertEdge(f.vertex[0], f.vertex[1], edges);
      InsertEdge(f.vertex[1], f.vertex[2], edges);
      InsertEdge(f.vertex[2], f.vertex[0], edges);
   }
   cout << "Edge count: " << edges.size() << '\n';
   std::vector<int> faceCounts(10, 0);
   for (EdgeMap::iterator i = edges.begin(); i != edges.end(); ++i)
   {
      ++faceCounts[i->second.faceCount];
   }
   for (int i = 1; i < 6; ++i)
   {
      cout << "Edges with " << i << " faces: " << faceCounts[i] << '\n';
   }
}

bool CompareVerts(const Vec3& a, const Vec3& b)
{
   return std::lexicographical_compare(&a.x, &a.x+3, &b.x, &b.x+3);
}

bool MeshAlgorithm::CheckValidMesh(const PolyMesh& aMesh)
{
   bool valid = true;
   size_t vertCount = aMesh.verts.size();
   const PolyhedronList& polys = *aMesh.polys;
   for (size_t i = 0; i < polys.size(); ++i)
   {
      const Polyhedron& p = polys[i];
      for (int j = 0; j < p.arity; ++j)
      {
         int v = p.vertex[j];
         if (!(v >= 0 && v < (int)vertCount))
         {
            valid = false;
         }
      }
   }
   std::vector<Vec3> verts;
   verts = *aMesh.verts;
   std::sort(verts.begin(), verts.end(), CompareVerts);
   std::vector<Vec3>::iterator iter = std::unique(verts.begin(), verts.end());
   int duplicates = verts.end() - iter;
   if (iter != verts.end())
   {
      valid = false;
   }
   return valid;
}

bool MeshAlgorithm::CheckValidMesh(const Mesh& aMesh)
{
   bool valid = true;
   size_t vertCount = aMesh.verts.size();
   for (size_t i = 0; i < aMesh.faces.size(); ++i)
   {
      const Face& f = aMesh.faces[i];
      if (!(f.vertex[0] >= 0 && f.vertex[0] < (int)vertCount
         && f.vertex[1] >= 0 && f.vertex[1] < (int)vertCount
         && f.vertex[2] >= 0 && f.vertex[2] < (int)vertCount))
      {
         valid = false;
      }
   }
   return valid;
}

void MeshAlgorithm::CheckValidFaces(Mesh& aMesh)
{
   
}

void MeshAlgorithm::SmoothMesh(Mesh&            aMesh,
                               VertexFaceLinks& aLinks,
                               float            aSmoothFactor)
{
   std::vector<Vec3> averages(aMesh.verts.size(), Vec3());
   std::vector<float>  averageCount(averages.size(), 0.0f);
   int starVertices[128];
   for (size_t i = 0; i < aMesh.verts.size(); ++i)
   {
      int vertCount = GetVertexStar(i, aMesh, aLinks, starVertices);
      for (int v = 0; v < vertCount; ++v)
      {
         int vertexIndex = starVertices[v];
         //float dist = (aMesh.verts[vertexIndex] - aMesh.verts[i]).length();
         float dist = 1.0f;
         if (dist > 0)
         {
            float distFactor = 1.0f / dist;
            averages[vertexIndex] += aMesh.verts[i] * distFactor;
            //++averageCount[vertexIndex];
            averageCount[vertexIndex] += distFactor;
         }
      }
   }
   for (size_t i = 0; i < averages.size(); ++i)
   {
      if (averageCount[i] > 0.0f)
      {
         Vec3 avg = averages[i] / averageCount[i];
         aMesh.verts[i] = interpolate(aSmoothFactor, aMesh.verts[i], avg);
      }
   }
}

void MeshAlgorithm::ComputeMeshParameters(const VertexList& aVerts,
                                          Box3&             aBounds)
{
   if (aVerts.empty())
   {
      aBounds = Box3();
   }
   else
   {
      float maxFloat = std::numeric_limits<float>::max();
      float minFloat = std::numeric_limits<float>::min();
      aBounds.Min() = Vec3(maxFloat, maxFloat, maxFloat);
      aBounds.Max() = Vec3(minFloat, minFloat, minFloat);
      for (size_t i = 0; i < aVerts.size(); ++i)
      {
         const Vec3& v = aVerts[i];
         for (int j = 0; j < 3; ++j)
         {
            aBounds.Min()[j] = std::min(v[j], aBounds.Min()[j]);
            aBounds.Max()[j] = std::max(v[j], aBounds.Max()[j]);
         }
      }
   }
}

void MeshAlgorithm::ComputeNonManifoldEdges(const Mesh&                   aMesh,
                                            std::vector<NonManifoldEdge>& aEdges)
{
   EdgeMap edges;
   for (size_t i = 0; i < aMesh.faces.size(); ++i)
   {
      const Face& f = aMesh.faces[i];
      InsertEdge(f.vertex[0], f.vertex[1], edges);
      InsertEdge(f.vertex[1], f.vertex[2], edges);
      InsertEdge(f.vertex[2], f.vertex[0], edges);
   }
   for (EdgeMap::iterator i = edges.begin(); i != edges.end(); ++i)
   {
      if (i->second.faceCount != 2)
      {
         NonManifoldEdge e;
         e.face_count = i->second.faceCount;
         e.verts[0] = i->first.first;
         e.verts[1] = i->first.second;
         aEdges.push_back(e);
      }
   }
}

void MeshAlgorithm::ComputeMaterialNonManifoldEdges(const Mesh&                   aMesh,
                                                    std::vector<NonManifoldEdge>& aEdges)
{
   MaterialEdgeMap edges;
   for (size_t i = 0; i < aMesh.faces.size(); ++i)
   {
      const Face& f = aMesh.faces[i];
      InsertEdge(f.vertex[0], f.vertex[1], f.materials[0], edges);
      InsertEdge(f.vertex[1], f.vertex[2], f.materials[0], edges);
      InsertEdge(f.vertex[2], f.vertex[0], f.materials[0], edges);
      InsertEdge(f.vertex[0], f.vertex[1], f.materials[1], edges);
      InsertEdge(f.vertex[1], f.vertex[2], f.materials[1], edges);
      InsertEdge(f.vertex[2], f.vertex[0], f.materials[1], edges);
   }
   for (MaterialEdgeMap::iterator i = edges.begin(); i != edges.end(); ++i)
   {
      if (i->second.faceCount != 2
         && (i->first.second == 3 || i->first.second == 4))
      {
         NonManifoldEdge e;
         e.face_count = i->second.faceCount;
         e.verts[0] = i->first.first.first;
         e.verts[1] = i->first.first.second;
         aEdges.push_back(e);
      }
   }
}

void MeshAlgorithm::AverageSmooth(Mesh& aMesh,
                                  int   aIterations)
{
   const int cMAX_IND = 128;
   size_t vertexCount = aMesh.verts.size();
   std::vector<Vec3> vertexMove(vertexCount);
   VertexFaceLinks links(cMAX_IND, aMesh);
   for (int iter = 0; iter < aIterations; ++iter)
   {
      for (size_t i = 0; i < vertexCount; ++i)
      {
         int indices[256];
         int n = GetVertexStar(i, aMesh, links, indices);
         int sumCount = 0;
         for (int v = 0; v < n; ++v)
         {
            int nVert = indices[v];
            const Vec3& vert = aMesh.verts[nVert];
            vertexMove[i] += vert;
            ++sumCount;
         }
         if (sumCount > 1)
         {
            vertexMove[i] /= (float)sumCount;
         }
         else
         {
            vertexMove[i] = aMesh.verts[i];
         }
      }

      for (size_t i = 0; i < vertexCount; ++i)
      {
         Vec3 avg = vertexMove[i];
         Vec3 start = aMesh.verts[i];
         Vec3 final = interpolate(.5, aMesh.verts[i], avg);
         aMesh.verts[i] = final;
         vertexMove[i] = Vec3();
      }
   }
}


float ComputeAngles(Mesh& aMesh,
                    Vec3& v1,
                    Vec3& v2,
                    Vec3& v3)
{
   Vec3 a(v1-v2), b(v3-v1);
   float cosA = dot(a, b) / (a.length() * b.length());
   float sinA = sqrt(1 - cosA * cosA);
   float cotA = -cosA / sinA;
   return cotA;
}

void MeshAlgorithm::CurvatureSmooth(Mesh& aMesh,
                                    int   aIterations)
{
   EdgeMesh em(aMesh);
   size_t vertexCount = aMesh.verts.size();
   //   std::vector<Vec3> vertexMove(vertexCount);
   std::vector<float> cots(em.GetEdges().size(), 0.0f);
   std::vector<Vec3>  vertexMove(vertexCount);
   std::vector<float> vertexScale(vertexCount, 0.0f);
   const EdgeMesh::FaceList& faces = em.GetFace();
   const EdgeMesh::EdgeList& edges = em.GetEdges();
   for (size_t i = 0; i < faces.size(); ++i)
   {
      const EdgeMeshData::EdgeFace& f = faces[i];
      Vec3& v1 = aMesh.verts[f.face.vertex[0]];
      Vec3& v2 = aMesh.verts[f.face.vertex[1]];
      Vec3& v3 = aMesh.verts[f.face.vertex[2]];
      cots[f.edges[2]] += ComputeAngles(aMesh, v1, v2, v3);
      cots[f.edges[0]] += ComputeAngles(aMesh, v2, v3, v1);
      cots[f.edges[1]] += ComputeAngles(aMesh, v3, v1, v2);
   }
   for (size_t i = 0; i < edges.size(); ++i)
   {
      const EdgeMeshData::EdgeEdge& e = edges[i];
      Vec3& v1 = aMesh.verts[e.verts[0]];
      Vec3& v2 = aMesh.verts[e.verts[1]];
      vertexMove[e.verts[0]] += cots[i] * v2;
      vertexMove[e.verts[1]] += cots[i] * v1;
      vertexScale[e.verts[0]] += cots[i];
      vertexScale[e.verts[1]] += cots[i];
   }
   for (size_t i = 0; i < vertexCount; ++i)
   {
      Vec3 avg = vertexMove[i] / vertexScale[i];
      Vec3 start = aMesh.verts[i];
      Vec3 final = interpolate(.5, aMesh.verts[i], avg);
      aMesh.verts[i] = final;
      vertexMove[i] = Vec3();
   }
}

void MeshAlgorithm::DistanceSmooth(Mesh& aMesh,
                                   int   aIterations)
{
   const int cMAX_IND = 128;
   size_t vertexCount = aMesh.verts.size();
   std::vector<Vec3> vertexMove(vertexCount);
   VertexFaceLinks links(cMAX_IND, aMesh);
   for (int iter = 0; iter < aIterations; ++iter)
   {
      for (size_t i = 0; i < vertexCount; ++i)
      {
         int indices[256];
         int n = GetVertexStar(i, aMesh, links, indices);
         float sumDist(0);
         const Vec3& cVert = aMesh.verts[i];
         for (int v = 0; v < n; ++v)
         {
            int nVert = indices[v];
            const Vec3& vert = aMesh.verts[nVert];
            float dist = (cVert - vert).length();
            //vertexMove[i] += vert;
            float mul = 1 / dist;
            vertexMove[i] += vert * mul;
            sumDist += mul;
         }
         if (sumDist > 0)
         {
            vertexMove[i] /= sumDist;
         }
         else
         {
            vertexMove[i] = aMesh.verts[i];
         }
      }

      for (size_t i = 0; i < vertexCount; ++i)
      {
         Vec3 avg = vertexMove[i];
         Vec3 start = aMesh.verts[i];
         Vec3 final = interpolate(.5, aMesh.verts[i], avg);
         aMesh.verts[i] = final;
         vertexMove[i] = Vec3();
      }
   }
}


void MeshAlgorithm::FlipMaterials(Mesh& aMesh)
{
   for (size_t i = 0; i < aMesh.faces.size(); ++i)
   {
      std::swap(aMesh.faces[i].materials[0], aMesh.faces[i].materials[1]);
   }
}

void MeshAlgorithm::RefineMesh(EdgeMesh& aMesh,
                               float     aLength)
{
   typedef std::pair<float, int> EdgeLen;
   std::priority_queue<EdgeLen> edgeQueue;
   float minLen = aLength * aLength;
   const std::vector<Vec3>& verts = aMesh.GetVerts();
   const EdgeMesh::EdgeList& edges = aMesh.GetEdges();
   for (size_t i = 0; i < edges.size(); ++i)
   {
      const EdgeMeshData::EdgeEdge& e = edges[i];
      float l2 = (verts[e.verts[0]] - verts[e.verts[1]]).length_squared();
      if (l2 > minLen)
      {
         edgeQueue.push(EdgeLen(l2, i));
      }
   }
   while (! edgeQueue.empty())
   {
      EdgeLen el = edgeQueue.top();
      edgeQueue.pop();
      std::vector<int> newEdges;
      aMesh.SplitEdge(el.second, newEdges);
      //for (size_t i = 0; i < newEdges.size(); ++i)
      //{
      //   int newId = newEdges[i];
      //   const EdgeMeshData::EdgeEdge& e = edges[newId];
      //   float l2 = (verts[e.verts[0]] - verts[e.verts[1]]).length_squared();
      //   if (l2 > minLen)
      //   {
      //      edgeQueue.push(EdgeLen(l2, newId));
      //   }
      //}
   }
}
void MeshAlgorithm::RefineMeshRatio(EdgeMesh& aMesh,
                                    float     aMaxRatio)
{
   const EdgeMesh::FaceList& faces = aMesh.GetFace();
   const std::vector<Vec3>& verts = aMesh.GetVerts();
   const EdgeMesh::EdgeList& edges = aMesh.GetEdges();

   std::vector<float>   edgeLengths(edges.size());
   std::vector<float>   edgeRatios(edges.size(), 0.0f);
   
   for (size_t i = 0; i < edges.size(); ++i)
   {
      const EdgeMeshData::EdgeEdge& e = edges[i];
      float l = (verts[e.verts[0]] - verts[e.verts[1]]).length();
      edgeLengths[i] = l;
   }

   for (size_t i = 0; i < faces.size(); ++i)
   {
      const EdgeMeshData::EdgeFace& f = faces[i];
      edgeRatios[f.edges[0]] = std::max(edgeRatios[f.edges[0]], 
         edgeLengths[f.edges[0]] / std::min(edgeLengths[f.edges[1]], edgeLengths[f.edges[2]]));
      edgeRatios[f.edges[1]] = std::max(edgeRatios[f.edges[1]], 
         edgeLengths[f.edges[1]] / std::min(edgeLengths[f.edges[0]], edgeLengths[f.edges[2]]));
      edgeRatios[f.edges[2]] = std::max(edgeRatios[f.edges[2]], 
         edgeLengths[f.edges[2]] / std::min(edgeLengths[f.edges[0]], edgeLengths[f.edges[1]]));
   }

   typedef std::pair<float, int> EdgeRatio;
   std::vector<EdgeRatio> ratios;
   for (size_t i = 0; i < edges.size(); ++i)
   {
      if (edgeRatios[i] > aMaxRatio)
      {
         ratios.push_back(EdgeRatio(edgeRatios[i], i));
      }
   }

   std::sort(ratios.begin(), ratios.end(), std::greater<EdgeRatio>());
   
   size_t splitCount = ratios.size();
   
   for (size_t i = 0; i < splitCount; ++i)
   {
      std::vector<int> newEdges;
      int edge = ratios[i].second;
      aMesh.SplitEdge(edge, newEdges);
   }
}

class SplitCriteria
{
public:
   SplitCriteria(Mesh& aMesh) : mMeshPtr(&aMesh) { }
   bool operator()(Face& aFace)
   {
      return mMeshPtr->verts[aFace.vertex[0]].x > 0;
   }
   Mesh* mMeshPtr;
};


void MeshAlgorithm::SplitMesh(Mesh& aMesh)
{
   SplitCriteria crit(aMesh);
   RemoveFaces(aMesh, crit);
}

void MeshAlgorithm::OrientMaterial(Mesh& aMesh,
                                   int   aMaterial)
{
   for (size_t i = 0; i < aMesh.faces.size(); ++i)
   {
      Face& f = aMesh.faces[i];
      if (f.materials[1] == aMaterial)
      {
         std::swap(f.materials[0], f.materials[1]);
         std::swap(f.vertex[0], f.vertex[2]);
      }
   }
}
typedef std::pair<int, pair<int, int> > SFace;
void ArrangeSFace(SFace& f)
{
   if (f.second.first > f.second.second)
   {
      swap(f.second.first, f.second.second);
   }
   if (f.first > f.second.first)
   {
      swap(f.first, f.second.first);
   }
   if (f.second.first > f.second.second)
   {
      swap(f.second.first, f.second.second);
   }
}
void MeshAlgorithm::FindDuplicateFaces(const Mesh&        aMesh,
                                       std::vector<Face>& aFaces)
{
   std::map<SFace, std::vector<Face> > badFaces;
   std::set<SFace> faces;
   for (size_t i = 0; i < aMesh.faces.size(); ++i)
   {
      const Face& f = aMesh.faces[i];
      SFace sf;
      sf.first = f.vertex[0];
      sf.second.first = f.vertex[1];
      sf.second.second = f.vertex[2];
      ArrangeSFace(sf);
      if (faces.find(sf) == faces.end())
      {
         faces.insert(sf);
      }
      else
      {
         aFaces.push_back(f);
         badFaces[sf];
      }
   }

   for (size_t i = 0; i < aMesh.faces.size(); ++i)
   {
      const Face& f = aMesh.faces[i];
      SFace sf;
      sf.first = f.vertex[0];
      sf.second.first = f.vertex[1];
      sf.second.second = f.vertex[2];
      ArrangeSFace(sf);
      if (badFaces.find(sf) != badFaces.end())
      {
         badFaces[sf].push_back(f);
      }
   }

   for (std::map<SFace, std::vector<Face> >::iterator iter = badFaces.begin(); iter != badFaces.end(); ++iter)
   {
      cout << "Dup faces: \n";
      for (size_t i = 0; i < iter->second.size(); ++i)
      {
         cout << "   " << i << ':';
         Face& f = iter->second[i];
         cout << " (" << f.vertex[0] << ',' << f.vertex[1] << ',' << f.vertex[2] << ')';
         cout << " = " << f.materials[0] << ' ' << f.materials[1] << '\n';
      }
   }
}

void MeshAlgorithm::CleanBadFaces(Mesh& aMesh)
{
   std::map<SFace, std::vector<int> > badFaces;
   std::set<SFace> faces;
   for (size_t i = 0; i < aMesh.faces.size(); ++i)
   {
      const Face& f = aMesh.faces[i];
      SFace sf;
      sf.first = f.vertex[0];
      sf.second.first = f.vertex[1];
      sf.second.second = f.vertex[2];
      ArrangeSFace(sf);
      if (faces.find(sf) == faces.end())
      {
         faces.insert(sf);
      }
      else
      {
         badFaces[sf];
      }
   }

   for (size_t i = 0; i < aMesh.faces.size(); ++i)
   {
      const Face& f = aMesh.faces[i];
      SFace sf;
      sf.first = f.vertex[0];
      sf.second.first = f.vertex[1];
      sf.second.second = f.vertex[2];
      ArrangeSFace(sf);
      if (badFaces.find(sf) != badFaces.end())
      {
         badFaces[sf].push_back(i);
      }
   }

   std::vector<int> facesToRemove;
   for (std::map<SFace, std::vector<int> >::iterator iter = badFaces.begin(); iter != badFaces.end(); ++iter)
   {
      cout << "Dup faces: \n";
      std::vector<int>& badFaces = iter->second;
      std::set<int> mats;
      for (size_t i = 0; i < badFaces.size(); ++i)
      {
         Face& f = aMesh.faces[badFaces[i]];
         mats.insert(f.materials[0]);
         mats.insert(f.materials[1]);
      }
      mats.erase(-1);
      if (mats.size() == 1)
      {
         mats.insert(-1);
      }
      if (mats.size() == 2)
      {
         Face& newFace = aMesh.faces[badFaces.front()];
         if (newFace.materials[0] != -1)
         {
            mats.erase(newFace.materials[0]);
            newFace.materials[1] = *mats.begin();
         }
         else if (newFace.materials[1] != -1)
         {
            mats.erase(newFace.materials[1]);
            newFace.materials[0] = *mats.begin();
         }
         facesToRemove.insert(facesToRemove.end(), badFaces.begin() + 1, badFaces.end());
      }
   }
   std::sort(facesToRemove.begin(), facesToRemove.end(), std::greater<int>());
   facesToRemove.erase(std::unique(facesToRemove.begin(), facesToRemove.end()), facesToRemove.end());
   if (! facesToRemove.empty())
   {
      cout << "Removing " << facesToRemove.size() << '\n';
   }
   //cout << "Erasing: ";
   for (size_t i = 0; i < facesToRemove.size(); ++i)
   {
      //cout << facesToRemove[i] << ", ";
      aMesh.faces->erase(aMesh.faces.begin() + facesToRemove[i]);
   }
   //FaceSelectionCriteria crit(facesToRemove);
   //RemoveFaces(aMesh, crit);
}


bool MeshAlgorithm::RayTriangleIntersect(Vec3 ray_origin, Vec3 ray_direction,
                                         Vec3 vert0, Vec3 vert1, Vec3 vert2,
                                         float& t, float& u, float& v)
{
   t = 0; u = 0; v = 0;

   Vec3 edge1 = vert1 - vert0;
   Vec3 edge2 = vert2 - vert0;

   Vec3 tvec, pvec, qvec;
   float det, inv_det;

   pvec = cross(ray_direction, edge2);

   det = dot(edge1, pvec);

   if (det > -0.00001f)
      return false;

   inv_det = 1.0f / det;

   tvec = ray_origin - vert0;

   u = dot(tvec, pvec) * inv_det;
   if (u < -0.001f || u > 1.001f)
      return false;

   qvec = cross(tvec, edge1);

   v = dot(ray_direction, qvec) * inv_det;
   if (v < -0.001f || u + v > 1.001f)
      return false;

   t = dot(edge2, qvec) * inv_det;

   if (t <= 0)
      return false;

   return true;
}

void MeshAlgorithm::FindRayHits(const Mesh& aMesh,
                                Vec3        aOrigin,
                                Vec3        aDirection,
                                std::vector<std::pair<float, int> >&     aHitFaces)
{
   for (size_t i = 0; i < aMesh.faces.size(); ++i)
   {
      float t,u,v;
      const Face& f = aMesh.faces[i];
      if (
         RayTriangleIntersect(aOrigin, aDirection, 
         aMesh.verts[f.vertex[0]],
         aMesh.verts[f.vertex[1]],
         aMesh.verts[f.vertex[2]], t,u,v)
         ||
         RayTriangleIntersect(aOrigin, aDirection, 
         aMesh.verts[f.vertex[2]],
         aMesh.verts[f.vertex[1]],
         aMesh.verts[f.vertex[0]], t,u,v)
         )
      {
         aHitFaces.push_back(pair<float, int>(t, i));
      }
   }
   std::sort(aHitFaces.begin(), aHitFaces.end());
}

void MeshAlgorithm::FindRayHits(const PolyMesh& aMesh,
                                Vec3              aOrigin,
                                Vec3              aDirection,
                                std::vector<std::pair<float, int> >& aHitPolys)
{
   for (size_t i = 0; i < aMesh.polys.size(); ++i)
   {
      float t,u,v;
      const Polyhedron& p = aMesh.polys[i];
      int faces = p.GetFaceCount();
      for (int j = 0; j < faces; ++j)
      {
         Polyhedron::FaceVerts f;
         p.GetFace(j, f);
         if (
            RayTriangleIntersect(aOrigin, aDirection, 
            aMesh.verts[f[0]],
            aMesh.verts[f[1]],
            aMesh.verts[f[2]], t,u,v)
            ||
            RayTriangleIntersect(aOrigin, aDirection, 
            aMesh.verts[f[2]],
            aMesh.verts[f[1]],
            aMesh.verts[f[0]], t,u,v)
            )
         {
            aHitPolys.push_back(pair<float, int>(t, i));
         }
      }
      
   }
   std::sort(aHitPolys.begin(), aHitPolys.end());
}

void SplitGroup(const Mesh&       aMesh,
                std::vector<int>& aGroup,
                std::vector<int>& aLessGroup,
                std::vector<int>& aMoreGroup)
{
   Vec3 minPos(1.0E10f, 1.0e10f, 1.0e10f), maxPos(-1.0E10f, -1.0E10f, -1.0E10f), avg;
   for (size_t i = 0; i < aGroup.size(); ++i)
   {
      int fIndex = aGroup[i];
      const Face& f = aMesh.faces[fIndex];
      const Vec3& a = aMesh.verts[f.vertex[0]];
      const Vec3& b = aMesh.verts[f.vertex[1]];
      const Vec3& c = aMesh.verts[f.vertex[2]];
      
      minPos = Vec3(std::min(minPos.x, a.x), std::min(minPos.y, a.y), std::min(minPos.z, a.z));
      minPos = Vec3(std::min(minPos.x, b.x), std::min(minPos.y, b.y), std::min(minPos.z, b.z));
      minPos = Vec3(std::min(minPos.x, c.x), std::min(minPos.y, c.y), std::min(minPos.z, c.z));

      maxPos = Vec3(std::max(maxPos.x, a.x), std::max(maxPos.y, a.y), std::max(maxPos.z, a.z));
      maxPos = Vec3(std::max(maxPos.x, b.x), std::max(maxPos.y, b.y), std::max(maxPos.z, b.z));
      maxPos = Vec3(std::max(maxPos.x, c.x), std::max(maxPos.y, c.y), std::max(maxPos.z, c.z));

      avg += a;
      avg += b;
      avg += c;
   }
   avg /= (float)(aGroup.size() * 3);
   int splitPlane = 0;
   Vec3 range = maxPos - minPos;
   if (range.y > range.x && range.y > range.z)
   {
      splitPlane = 1;
   }
   if (range.z > range.x && range.z > range.y)
   {
      splitPlane = 2;
   }
   float splitLocation = avg[splitPlane];
   for (size_t i = 0; i < aGroup.size(); ++i)
   {
      int fIndex = aGroup[i];
      const Face& f = aMesh.faces[fIndex];
      bool isLess = aMesh.verts[f.vertex[0]][splitPlane] <= splitLocation
         || aMesh.verts[f.vertex[1]][splitPlane] <= splitLocation
         || aMesh.verts[f.vertex[2]][splitPlane] <= splitLocation;
      bool isMore = aMesh.verts[f.vertex[0]][splitPlane] > splitLocation
         || aMesh.verts[f.vertex[1]][splitPlane] > splitLocation
         || aMesh.verts[f.vertex[2]][splitPlane] > splitLocation;
      if (isLess)
      {
         aLessGroup.push_back(fIndex);
      }
      if (isMore)
      {
         aMoreGroup.push_back(fIndex);
      }
   }
}

void MeshAlgorithm::DivideArea(const Mesh&                   aMesh,
                               std::list<std::vector<int> >& aFaceGroups)
{
   const int cMIN_GROUP_SIZE = 10;
   const float cMIN_SPLIT_FRACTION = 1/3.0f;
   typedef std::vector<int> FaceGroup;
   std::list<FaceGroup>& groups(aFaceGroups);
   
   groups.push_back(FaceGroup());
   FaceGroup& rootGroup = groups.back();
   for (size_t i = 0; i < aMesh.faces.size(); ++i)
   {
      rootGroup.push_back(i);
   }

   std::list<FaceGroup>::iterator groupIter = groups.begin();
   for (;groupIter != groups.end();)
   {
      FaceGroup& grp = *groupIter;
      FaceGroup lGrp, mGrp;
      if (grp.size() > cMIN_GROUP_SIZE)
      {
         SplitGroup(aMesh, grp, lGrp, mGrp);
         float frac = (lGrp.size() + mGrp.size() - grp.size()) / float(grp.size());
         if (frac < cMIN_SPLIT_FRACTION)
         {
            groups.push_back(FaceGroup());
            std::swap(groups.back(), lGrp);
            groups.push_back(FaceGroup());
            std::swap(groups.back(), mGrp);
            groupIter = groups.erase(groupIter);
         }
         else
         {
            ++groupIter;
         }
      }
      else
      {
         ++groupIter;
      }
   }
}

int WhichSide(Vec3 p1, Vec3 p2, Vec3 p3, Vec3 D, Vec3 P)
{
   int pos(0), neg(0);
   float t = D.dot(p1-P);
   if (t > 0) pos++; else if (t < 0) neg++;
   if (pos&&neg) return 0;
   t = D.dot(p2-P);
   if (t > 0) pos++; else if (t < 0) neg++;
   if (pos&&neg) return 0;
   t = D.dot(p3-P);
   if (t > 0) pos++; else if (t < 0) neg++;
   if (pos&&neg) return 0;
   return pos ? 1 : -1;
}

int tri_tri_intersect(float V0[3],float V1[3],float V2[3],
                      float U0[3],float U1[3],float U2[3]);
bool TriangleIntersect(Vec3 a1, Vec3 a2, Vec3 a3, Vec3 b1, Vec3 b2, Vec3 b3)
{
   //Vec3 x = cross(b3-b2, b1-b2);
   //return WhichSide(a1, a2, a3, x, b2) == 0;
   return 0!=tri_tri_intersect(&a1.x, &a2.x, &a3.x, &b1.x, &b2.x, &b3.x);
}

void MeshAlgorithm::FindIntersectingFaces(const Mesh&       aMesh,
                                          std::vector<int>& aFaces)
{
   typedef std::vector<int> FaceGroup;
   std::list<FaceGroup > faceGroups;
   DivideArea(aMesh, faceGroups);
   for (std::list<FaceGroup>::iterator i = faceGroups.begin(); i != faceGroups.end(); ++i)
   {
      FaceGroup& grp = *i;
      for (size_t x = 0; x < grp.size(); ++x)
      {
         const Face& f1 = aMesh.faces[grp[x]];
         for (size_t y = x+1; y < grp.size(); ++y)
         {
            const Face& f2 = aMesh.faces[grp[y]];
            bool hit = false;
            if (
             f2.vertex[0] != f1.vertex[0]
            && f2.vertex[0] != f1.vertex[1]
            && f2.vertex[0] != f1.vertex[2]
            && f2.vertex[1] != f1.vertex[0]
            && f2.vertex[1] != f1.vertex[1]
            && f2.vertex[1] != f1.vertex[2]
            && f2.vertex[2] != f1.vertex[0]
            && f2.vertex[2] != f1.vertex[1]
            && f2.vertex[2] != f1.vertex[2]
            )
            {
               hit = TriangleIntersect(
               aMesh.verts[f2.vertex[0]],
               aMesh.verts[f2.vertex[1]],
               aMesh.verts[f2.vertex[2]],
               aMesh.verts[f1.vertex[0]],
               aMesh.verts[f1.vertex[1]], 
               aMesh.verts[f1.vertex[2]]);
            }
            if (hit)
            {
               aFaces.push_back(grp[x]);
               aFaces.push_back(grp[y]);
            }
         }
      }
   }
   std::sort(aFaces.begin(), aFaces.end());
   aFaces.erase(std::unique(aFaces.begin(), aFaces.end()), aFaces.end());
}

void MeshAlgorithm::BoxSelectVertices(std::vector<Vec3>& aVertexList,
                                      Mat4               aTransform,
                                      Vec2               aMin,
                                      Vec2               aMax,
                                      std::vector<int>&  aInside)
{
   for (size_t i = 0; i < aVertexList.size(); ++i)
   {
      Vec3 v = aVertexList[i] * aTransform;
      if (v.x >= aMin.x && v.x <= aMax.x && v.y>=aMin.y && v.y <= aMax.y)
      {
         aInside.push_back(i);
      }
   }
}

int AddMidpoint(PolyMesh& aMesh, std::map<Edge, int>& aMidpoints, int i1, int i2)
{
   int vert;
   Edge e = OrderEdge(i1, i2);
   std::map<Edge, int>::iterator iter = aMidpoints.find(e);
   if (iter == aMidpoints.end())
   {
      vert = aMidpoints[e] = aMesh.verts.size();
      Vec3 mid = (aMesh.verts[i1] + aMesh.verts[i2]) / 2.0f;
      aMesh.verts.push_back(mid);
   }
   else
   {
      vert = iter->second;
   }
   return vert;
}
//void MeshAlgorithm::Subdivide(PolyMesh& aMesh)
//{
//   std::map<Edge, int> midpoints;
//   for (size_t i = 0; i < aMesh.tetra.size(); ++i)
//   {
//      Tetra& t = aMesh.tetra[i];
//      int v0 = t.vertex[0];
//      int v1 = t.vertex[1];
//      int v2 = t.vertex[2];
//      int v3 = t.vertex[3];
//      int m01 = AddMidpoint(aMesh, midpoints, t.vertex[0], t.vertex[1]);
//      int m02 = AddMidpoint(aMesh, midpoints, t.vertex[0], t.vertex[2]);
//      int m03 = AddMidpoint(aMesh, midpoints, t.vertex[0], t.vertex[3]);
//      int m12 = AddMidpoint(aMesh, midpoints, t.vertex[1], t.vertex[2]);
//      int m13 = AddMidpoint(aMesh, midpoints, t.vertex[1], t.vertex[3]);
//      int m23 = AddMidpoint(aMesh, midpoints, t.vertex[2], t.vertex[3]);
//
//      int tets[4][4] = {
//               { v0, m01, m02, m03 },
//               { v1, m01, m12, m13 },
//               { v2, m02, m12, m23 },
//               { v3, m03, m13, m23 }
//      };
//      for (int ti = 0; ti < 4; ++ti)
//      {
//         Tetra nt;
//         nt.material = t.material;
//         std::copy(tets[ti], tets[ti+1], nt.vertex.begin());
//         aMesh.tetra.push_back(nt);
//      }
//      Octa oct;
//      oct.material = t.material;
//      int octVerts[6] = { m01, m12, m02, m03, m13, m23 };
//      std::copy(octVerts, octVerts + 6, oct.vertex.begin());
//      aMesh.octa.push_back(oct);
//   }
//}

struct EdgeFaceData
{
   int               startVertex;
   int               nextVertex;
   std::vector<int>  tris;
};
class CycleFixer
{
   std::map<int, int> toNewVert;
   Mesh* mMesh;
   AdjacencyMap mAdjacencyMap;
   int NoDupVert(int aIndex)
   {
      toNewVert[aIndex] = aIndex;
      return aIndex;
   }
   int DupVert(int aIndex)
   {
      std::map<int,int>::iterator i = toNewVert.find(aIndex);
      if (i != toNewVert.end()) { return i->second; }
      else
      {
         int v = toNewVert[aIndex] = mMesh->verts.size();
         Vec3 dv = mMesh->verts[aIndex];
         //dv.z += -.05;
         mMesh->verts.push_back(dv);
         return v;
      }
   }
   //struct RegionData
   //{
   //   RegionInterface mRegion;
   //   EdgeList
   //};
public:
   
   struct LoopData
   {
      RegionInterface   mInterface;
      int               mInteriorTris;
      VertexLoop::ID    mLoopId;
      bool operator<(const LoopData& aRhs) const
      {
         if (mLoopId < aRhs.mLoopId) return true;
         if (aRhs.mLoopId < mLoopId) return false;
         return mInteriorTris < aRhs.mInteriorTris;
      }
   };
   bool FindCycles(PolyMesh&                     aMesh,
                   CreaseStructure&              aCreases,
                   std::vector<int>&             aSelVerts,
                   std::list<std::pair<RegionInterface, IndexList> >& aCycles)
   {
      bool fixed(false);
      std::set<int> selVerts(aSelVerts.begin(), aSelVerts.end());
      std::vector<int> selLoops;
      for (size_t i = 0; i < aCreases.mLoops.size(); ++i)
      {
         VertexLoop& loop = *aCreases.mLoops[i];
         for (size_t j = 0; j < loop.mVertices.size(); ++j)
         {
            if (selVerts.find(loop.mVertices[j]) != selVerts.end())
            {
               int bestRegion = 0;
               size_t bestFaces = 999999999;
               for (size_t k = 0; k < loop.mCreases.size(); ++k)
               {
                  RegionCrease* c = loop.mCreases[k];
                  if (c->mFaces.size() < bestFaces)
                  {
                     bestRegion = k;
                     bestFaces = c->mFaces.size();
                  }
               }
               FixCycle(aMesh.GetCreaseMesh(), aCreases, i, loop.mCreases[bestRegion]->mInterface);
               return true;
            }
         }
      }
      return fixed;
   }

   bool FixCycle(Mesh&             aMesh,
                 CreaseStructure&  aCreases,
                 int               aLoopIndex,
                 RegionInterface   aInterface)
   {
      mMesh = &aMesh;
      MeshAlgorithm::MakeAdjacencyMap(*aMesh.faces, mAdjacencyMap);

      std::vector<EdgeFaceData> edges;
      std::vector<int> touchingTris;
      const VertexList& verts = *aMesh.verts;
      VertexLoop& cycleLoop = *aCreases.mLoops[aLoopIndex];
      
      // Find crease to be duplicated
      RegionCrease* dupCrease = 0;
      for (size_t i = 0; i < cycleLoop.mCreases.size(); ++i)
      {
         if (cycleLoop.mCreases[i]->mInterface == aInterface)
         {
            dupCrease = cycleLoop.mCreases[i];
            break;
         }
      }

      std::set<int> vertsOnLoop(cycleLoop.mVertices.begin(), cycleLoop.mVertices.end());
      FaceList& faces = *aMesh.faces;
      // Find triangles that share >= 2 vertices on this loop
      for (size_t i = 0; i < faces.size(); ++i)
      {
         const Face& f = faces[i];
         int vcount = (vertsOnLoop.find(f.vertex[0]) != vertsOnLoop.end()) ? 1 : 0;
         vcount += (vertsOnLoop.find(f.vertex[1]) != vertsOnLoop.end()) ? 1 : 0;
         if (vcount == 0)
         {
            vcount += (vertsOnLoop.find(f.vertex[2]) != vertsOnLoop.end()) ? 1 : 0;
         }
         if (vcount > 0)
         {
            touchingTris.push_back(i);
         }
      }

      // Collect edges and triangles that share the edge
      for (size_t i = 0; i < cycleLoop.mVertices.size(); ++i)
      {
         EdgeFaceData edge;
         int vert = edge.startVertex = cycleLoop.mVertices[i];
         int nextVert = edge.nextVertex = cycleLoop.mVertices[(i+1)%cycleLoop.mVertices.size()];
         for (size_t j = 0; j < touchingTris.size(); ++j)
         {
            Face& f = faces[touchingTris[j]];
            if (  (f.vertex[0] == vert && f.vertex[1] == nextVert)
               || (f.vertex[1] == vert && f.vertex[2] == nextVert)
               || (f.vertex[2] == vert && f.vertex[0] == nextVert)
               || (f.vertex[0] == nextVert && f.vertex[1] == vert)
               || (f.vertex[1] == nextVert && f.vertex[2] == vert)
               || (f.vertex[2] == nextVert && f.vertex[0] == vert)
               )
            {
               edge.tris.push_back(touchingTris[j]);
               //OrderFace(f, vert, nextVert);
            }
         }
         edges.push_back(edge);
      }

      int arity = edges.front().tris.size();

      // Make sure the loop has consistant arity
      {
         for (size_t i = 1; i < edges.size(); ++i)
         {
            if (edges[i].tris.size() != arity)
            {
               cout << "Failed fixing loop -- differing face count on edges\n";
               return false;
            }
         }
      }

      // Compute correct region ordering
      std::vector<RegionInterface> regionOrdering;
      {
         int firstRegion = 0;
         EdgeFaceData& firstEdge = edges.front();
         for (size_t i = 0; i < firstEdge.tris.size(); ++i)
         {
            Face& f = faces[firstEdge.tris[i]];
            OrderFace(f, firstEdge.startVertex, firstEdge.nextVertex);
         }
         for (size_t i = 0; i < firstEdge.tris.size(); ++i)
         {
            int ti = firstEdge.tris[i];
            const Face& f = faces[ti];
            if (region_equal(aInterface, f))
            {
               firstRegion = i;
               break;
            }
         }
         Face& firstFace = faces[firstEdge.tris[firstRegion]];
         Vec3 edgeV = (verts[firstEdge.nextVertex] - verts[firstEdge.startVertex]).normal();
         Vec3 midpt = (verts[firstEdge.startVertex] + verts[firstEdge.nextVertex]) / 2.0f;
         Vec3 toFirst = (verts[firstFace.vertex[2]] - midpt).normal();
         Vec3 toC = cross(toFirst, edgeV).normal();
         typedef std::pair<float, int> FaceAngle;
         std::vector<FaceAngle> faceAngles;
         for (size_t i = 0; i < firstEdge.tris.size(); ++i)
         {
            float angle(0.0f);
            if (i != firstRegion)
            {
               int ti = firstEdge.tris[i];
               const Face& f = faces[ti];
               Vec3 toFace = (verts[f.vertex[2]] - midpt).normal();
               angle = atan2(toFace.dot(toC), toFace.dot(toFirst));
            }
            if (angle < 0.0f) angle += Math::cTWO_PI;
            faceAngles.push_back(FaceAngle(angle, i));
         }
         std::sort(faceAngles.begin(), faceAngles.end());
         for (size_t i = 0; i < faceAngles.size(); ++i)
         {
            int fi = firstEdge.tris[faceAngles[i].second];
            regionOrdering.push_back(face_region(faces[fi]));
         }
      }

      {
         cout << "Regions: ";
         for (size_t i = 0; i < regionOrdering.size(); ++i)
         {
            RegionInterface interface = regionOrdering[i];
            cout << '(' << interface.first << ',' << interface.second << ')';
         }
         cout << '\n';
      }
      if (regionOrdering.size() <= 1)
      {
         return false;
      }

      // Order edge faces using region ordering
      {
         for (size_t i = 0; i < edges.size(); ++i)
         {
            EdgeFaceData& data = edges[i];
            std::vector<int> newTris(data.tris.size());
            for (size_t r = 0; r < regionOrdering.size(); ++r)
            {
               for (size_t j = 0; j < data.tris.size(); ++j)
               {
                  const Face& f = faces[data.tris[j]];
                  if (region_equal(regionOrdering[r], f))
                  {
                     newTris[r] = data.tris[j];
                     break;
                  }
               }
            }
            std::swap(newTris, data.tris);
         }
      }

      // Dont duplicate vertices on other loops in the duplicated region
      for (size_t i = 0; i < dupCrease->mLoops.size(); ++i)
      {
         VertexLoop* loopPtr = dupCrease->mLoops[i];
         if (loopPtr != &cycleLoop)
         {
            for (size_t j = 0; j < loopPtr->mVertices.size(); ++j)
            {
               NoDupVert(loopPtr->mVertices[j]);
            }
         }
      }

      // Figure out how to assign regions for new faces
      int innerRegion;
      int regionA;
      {
         Face& f = faces[edges.front().tris[1]];
         if (f.materials[0] == aInterface.first || f.materials[0] == aInterface.second)
         {
            innerRegion = f.materials[1];
            regionA = f.materials[0];
         }
         else 
         {
            innerRegion = f.materials[0];
            regionA = f.materials[1];
         }
      }
      cout << "Filling loop with material " << innerRegion << '\n';

      RegionInterface borderInterface = regionOrdering[1];
      for (size_t i = 0; i < touchingTris.size(); ++i)
      {
         Face& f = faces[touchingTris[i]];
         if (region_equal(borderInterface, f))
         {
            for (int v = 0; v < 3; ++v)
            {
               int vi = f.vertex[v];
               if (vertsOnLoop.find(vi) != vertsOnLoop.end())
               {
                  f.vertex[v] = DupVert(vi);
               }
            }
         }
      }
      
      // duplicate faces in region
      std::vector<int>& interiorFaces = dupCrease->mFaces;
      //std::vector<int> startfaces;
      //startfaces.push_back(edges.front().tris.front());
      //FillRegion(mAdjacencyMap, *mMesh->faces, startfaces, interiorFaces);
      cout << "Filling interior: " << interiorFaces.size() << " faces. \n";
      //int regionA = aInterface.first;
      //int innerRegion = newMaterial;
      //int regionC = aInterface.second;
      for (size_t i = 0; i < interiorFaces.size(); ++i)
      {
         Face f = faces[interiorFaces[i]];
         f.vertex[0] = DupVert(f.vertex[0]);
         f.vertex[1] = DupVert(f.vertex[1]);
         f.vertex[2] = DupVert(f.vertex[2]);
         ///!!! Selecting the correct materials ???
         //f.materials[1] = newMaterial;
         //faces[interiorFaces[i]].materials[0] = newMaterial;
         // duplicated face?
         if (f.materials[0] == regionA)
         {
            f.materials[1] = innerRegion;
            faces[interiorFaces[i]].materials[0] = innerRegion;
         }
         else 
         {
            f.materials[0] = innerRegion;
            faces[interiorFaces[i]].materials[1] = innerRegion;
         }
         faces.push_back(f);
      }
      return true;
   }

};
void MeshAlgorithm::FindCycles(PolyMesh&                     aMesh,
                               CreaseStructure&              aCreases,
                               std::vector<int>&             aSelVerts,
                std::list<std::pair<RegionInterface, IndexList> >& aCycles)
{
   CycleFixer f;
   f.FindCycles(aMesh, aCreases, aSelVerts, aCycles);
}


bool MeshAlgorithm::FixCycle(Mesh&            aMesh,
                             CreaseStructure& aCreases,
                             int              aLoopIndex)
{
   VertexLoop& loop = *aCreases.mLoops[aLoopIndex];
   int bestPatch=-1;
   size_t bestTris = 9000000;
   for (size_t i = 0; i < loop.mCreases.size(); ++i)
   {
      if (loop.mCreases[i]->mFaces.size() < bestTris)
      {
         bestTris = loop.mCreases[i]->mFaces.size();
         bestPatch = i;
      }
   }
   if (bestPatch != -1)
   {
      return FixCycle(aMesh, aCreases, aLoopIndex, loop.mCreases[bestPatch]->mInterface);
   }
   return false;
}

bool MeshAlgorithm::FixCycle(Mesh&             aMesh,
                             CreaseStructure&  aCreases,
                             int               aLoopIndex,
                             RegionInterface   aInterface)
{
   CycleFixer f;
   return f.FixCycle(aMesh, aCreases, aLoopIndex, aInterface);
}

bool MeshAlgorithm::DetatchInterface(Mesh&            aMesh,
                                     CreaseStructure& aCreases,
                                     int              aInterfaceIndex,
                                     int              aFillMaterial)
{
   return false;
}

void MeshAlgorithm::MakeAdjacencyMap(const FaceList&   aFaces,
                                     AdjacencyMap& aAdjacencyMap)
{
   DEBUG_LOG("MakeAdjacencyMap")
   for (size_t i = 0; i < aFaces.size(); ++i)
   {
      const Face& f = aFaces[i];
      Edge e[] = { Edge(f.vertex[0], f.vertex[1]), Edge(f.vertex[1], f.vertex[2]), Edge(f.vertex[2], f.vertex[0]) };
      for (int j = 0; j < 3; ++j)
      {
         OrderEdge(e[j]);
         aAdjacencyMap[e[j]].push_back(i);
      }
   }
   DEBUG_LOG(".\n")
}

void MeshAlgorithm::FindFaces(const FaceList&         aFaces,
                              const std::vector<int>& aVerts,
                              std::vector<int>&       aSelFaces)
{
   std::set<int> verts(aVerts.begin(), aVerts.end());
   for (size_t i = 0; i < aFaces.size(); ++i)
   {
      bool fail=false;
      const Face& f = aFaces[i];
      for (size_t j = 0; j < 3; ++j)
      {
         if (verts.find(f.vertex[j]) == verts.end()) { fail = true; break; }
      }
      if (!fail) aSelFaces.push_back(i);
   }
}

void MeshAlgorithm::FillRegion(MeshAdjacency&          aAdjacency,
                               const std::vector<int>& aStartFaces,
                               std::vector<int>&       aRegionFaces)
{
   if (aStartFaces.empty()) return;
   Mesh& mesh = aAdjacency.GetMesh();
   FaceList& aFaces = *mesh.faces;
   std::vector<int> nextFaces(aStartFaces);
   int fi = nextFaces.front();
   const Face& rootFace = mesh.faces[fi];
   RegionInterface interface = face_region(rootFace);
   OrderEdge(interface);
   std::set<int> visitedFaces;
   while (! nextFaces.empty())
   {
      int faceIndex = nextFaces.back();
      nextFaces.pop_back();
      if (visitedFaces.find(faceIndex) != visitedFaces.end()) continue;
      aRegionFaces.push_back(faceIndex);
      visitedFaces.insert(faceIndex);
      //visitedFaces.insert(faceIndex);
      const Face& f = aFaces[faceIndex];
      Edge e[] = { Edge(f.vertex[0], f.vertex[1]), Edge(f.vertex[1], f.vertex[2]), Edge(f.vertex[2], f.vertex[0]) };
      for (int j = 0; j < 3; ++j)
      {
         int faces[64];
         size_t fCount = aAdjacency.GetEdgeFaces(e[j], faces);
         for (size_t k = 0; k < fCount; ++k)
         {
            int ind = faces[k];
            const Face& nextFace = aFaces[ind];
            if (ind != faceIndex && region_equal(RegionInterface(nextFace.materials[0], nextFace.materials[1]), interface))
            {
               nextFaces.push_back(ind);
            }
         }
      }
   }
}


void MeshAlgorithm::FillRegion(MeshAlgorithm::AdjacencyMap& aAdjacencyMap,
                               const FaceList&                    aFaces,
                               const std::vector<int>&            aStartFaces,
                               std::vector<int>&                  aRegionFaces)
{
   std::vector<int> nextFaces(aStartFaces);
   int fi = nextFaces.front();
   const Face& rootFace = aFaces[fi];
   RegionInterface interface(rootFace.materials[0], rootFace.materials[1]);
   OrderEdge(interface);
   std::set<int> visitedFaces;
   while (! nextFaces.empty())
   {
      int faceIndex = nextFaces.back();
      nextFaces.pop_back();
      if (visitedFaces.find(faceIndex) != visitedFaces.end()) continue;
      aRegionFaces.push_back(faceIndex);
      visitedFaces.insert(faceIndex);
      //visitedFaces.insert(faceIndex);
      const Face& f = aFaces[faceIndex];
      Edge e[] = { Edge(f.vertex[0], f.vertex[1]), Edge(f.vertex[1], f.vertex[2]), Edge(f.vertex[2], f.vertex[0]) };
      for (int j = 0; j < 3; ++j)
      {
         OrderEdge(e[j]);
         const AdjacentFaces& a_faces = aAdjacencyMap[e[j]];
         for (size_t k = 0; k < a_faces.size(); ++k)
         {
            int ind = a_faces[k];
            const Face& nextFace = aFaces[ind];
            if (ind != faceIndex && region_equal(RegionInterface(nextFace.materials[0], nextFace.materials[1]), interface))
            {
               nextFaces.push_back(ind);
            }
         }
      }
   }
}

//
//struct RegionInfo
//{
//   RegionInterface mRegion;
//   int             mFaceCount;
//   EdgeList        mEdges;
//};
//
//void MeshAlgorithm::FindCreaseRegions(const Mesh&              aMesh,
//                                      std::list<RegionCrease>& aCreases)
//{
//   // Collect region info in this temporary storage
//   std::vector<RegionInfo> regionInfos;
//   
//   // build adjacency map
//   AdjacencyMap edges;
//   const std::vector<Face>& faces = *aMesh.faces;
//   MakeAdjacencyMap(faces, edges);
//
//   std::set<int> unvisitedFaces;
//   unvisitedFaces.insert(boost::counting_iterator<int>(0), boost::counting_iterator<int>(faces.size()));
//   while (! unvisitedFaces.empty())
//   {
//      std::vector<Edge> containedEdges;
//      //std::set<int> visitedFaces;
//      //std::vector<Edge> nextEdges;
//      std::vector<int> nextFaces;
//      std::vector<int> interiorFaces;
//      int fi = *unvisitedFaces.begin();
//      const Face& rootFace = faces[fi];
//      RegionInterface interface(rootFace.materials[0], rootFace.materials[1]);
//      OrderEdge(interface);
//      nextFaces.push_back(fi);
//      while (! nextFaces.empty())
//      {
//         int faceIndex = nextFaces.back();
//         nextFaces.pop_back();
//         if (unvisitedFaces.find(faceIndex) == unvisitedFaces.end()) { continue; }
//         interiorFaces.push_back(faceIndex);
//         unvisitedFaces.erase(faceIndex);
//         //visitedFaces.insert(faceIndex);
//         const Face& f = faces[faceIndex];
//         Edge e[] = { Edge(f.vertex[0], f.vertex[1]), Edge(f.vertex[1], f.vertex[2]), Edge(f.vertex[2], f.vertex[0]) };
//         for (int j = 0; j < 3; ++j)
//         {
//            OrderEdge(e[j]);
//            containedEdges.push_back(e[j]);
//            AdjacentFaces& a_faces = edges[e[j]];
//            for (size_t k = 0; k < a_faces.size(); ++k)
//            {
//               int ind = a_faces[k];
//               const Face& nextFace = faces[ind];
//               if (ind != faceIndex && region_equal(RegionInterface(nextFace.materials[0], nextFace.materials[1]), interface))
//               {
//                  nextFaces.push_back(ind);
//               }
//            }
//         }
//      }
//      std::vector<Edge> creases;
//      std::sort(containedEdges.begin(), containedEdges.end());
//      size_t i(1);
//      for (; i < containedEdges.size();)
//      {
//         if (containedEdges[i-1] != containedEdges[i])
//         {
//            creases.push_back(containedEdges[i-1]);
//            ++i;
//         }
//         else
//         {
//            while (i < containedEdges.size() && containedEdges[i] == containedEdges[i-1])
//            {
//               ++i;
//            }
//            ++i;
//         }
//      }
//      if (i == containedEdges.size())
//      {
//         creases.push_back(containedEdges[i-1]);
//      }
//      RegionInfo info;
//      info.mEdges = creases;
//      info.mRegion = interface;
//      info.mFaceCount = interiorFaces.size();
//      regionInfos.push_back(info);
//      //RegionCreases& regionCreases = mRegionCreases[interface];
//      //regionCreases.mRegions = interface;
//      //regionCreases.mCreases.push_back(EdgeList());
//      //regionCreases.mCreases.back().swap(creases);
//      //regionCreases.mInteriorFaceCount = interiorFaces.size();
//   }
//
//   // Now move from temporary storage to final structure
//   for (size_t i = 0; i < regionInfos.size(); ++i)
//   {
//      aCreases.push_back(RegionCrease());
//      RegionCrease& crease = aCreases.back();
//      RegionInfo& info = regionInfos[i];
//      crease.mRegions = info.mRegion;
//      crease.mInteriorFaceCount = info.mFaceCount;
//      EdgeCycleList cycles;
//      FindCycles(info.mEdges, 1000000, cycles);
//      if (cycles.size() > 0)
//      {
//         //crease.mLoops.resize(cycles.size());
//         for (EdgeCycleList::iterator iter = cycles.begin(); iter != cycles.end(); ++iter)
//         {
//            EdgeCycle& cycle = *iter;
//            crease.mLoops.push_back(VertexLoop());
//            VertexLoop& loop = crease.mLoops.back();
//            std::swap(loop.vertices, cycle);
//            NormalizeLoop(loop);
//            loop.RecomputeHash();
//         }
//      }
//      else
//      {
//         aCreases.pop_back();
//      }
//   }
//}

void MeshAlgorithm::GetTetPairs(PolyMesh&             aMesh,
                                std::vector<TetPair>& aTetPairs)
{
   typedef boost::array<int, 3> FaceInd;
   typedef std::map<FaceInd, int> TetFaceMap;
   TetFaceMap tetFaces;
   for (size_t i = 0; i < aMesh.polys.size(); ++i)
   {
      Polyhedron& poly = aMesh.polys[i];
      FaceInd verts;
      for (int j = 0; j < poly.GetFaceCount(); ++j)
      {
         poly.GetFace(j, verts);
         int tipIndex = poly.GetTip(j);
         OrderFace(verts);
         TetFaceMap::iterator iter = tetFaces.find(verts);
         if (iter != tetFaces.end())
         {
            TetPair p;
            p.tetTips[0] = iter->second;
            p.tetTips[1] = tipIndex;
            p.verts = verts;
            aTetPairs.push_back(p);
         }
         else
         {
            tetFaces[verts] = tipIndex;
         }
      }
   }
}

void MeshAlgorithm::ExtractMatPoints(Mesh&                  aMesh,
                                     std::vector<MatPoint>& aPoints)
{
   aPoints.assign(aMesh.verts.size(), MatPoint());
   for (size_t i = 0; i < aMesh.verts.size(); ++i)
   {
      aPoints[i].point = aMesh.verts[i];
   }
   for (size_t i = 0; i < aMesh.faces.size(); ++i)
   {
      Face& f = aMesh.faces[i];
      for (size_t j = 0; j < 3; ++j)
      {
         int vi = f.vertex[j];
         for (size_t m = 0; m < 2; ++m)
         {
            int mat = f.materials[m];
            tiny_vector<int, 3>& mats = aPoints[vi].materials;
            if (std::find(mats.begin(), mats.end(), mat) == mats.end())
            {
               mats.push_back(mat);
            }
         }
      }
   }
   for (size_t i = 0; i < aPoints.size(); ++i)
   {
      MatPoint& pt = aPoints[i];
      std::sort(pt.materials.begin(), pt.materials.end());
   }
}

void MeshAlgorithm::GetBorderEdges(const Mesh&      aMesh,
                                   const IndexList& aFaces,
                                   EdgeList&        aBorderEdges)
{
   std::map<Edge, int> edgeUses;
   for (size_t i = 0; i < aFaces.size(); ++i)
   {
      const Face& f = aMesh.faces[aFaces[i]];
      for (int ei = 0; ei < 3; ++ei)
      {
         ++edgeUses[f.GetEdge(ei)];
      }
   }
   for (std::map<Edge, int>::iterator i = edgeUses.begin(); i != edgeUses.end(); ++i)
   {
      if (i->second == 1)
      {
         aBorderEdges.push_back(i->first);
      }
   }
}

void MeshAlgorithm::GetEdgeVertices(const EdgeList& aEdges,
                                    IndexList&      aVertices)
{
   for (size_t i = 0; i < aEdges.size(); ++i)
   {
      aVertices.push_back(aEdges[i].first);
      aVertices.push_back(aEdges[i].second);
   }
   make_unique(aVertices);
}

void MeshAlgorithm::DetatchInterface4(Mesh&            aMesh,
                                      MeshAdjacency&   aAdjacency,
                                      const IndexList& aFaces,
                                      const EdgeList&  aEdges,
                                      int              aFillMaterial)
{
   int cutRegion(-999), otherRegion;
   if (! aFaces.empty())
   {
      Face& f = aMesh.faces[aFaces[0]];
      cutRegion = f.materials[0];
      otherRegion = f.materials[1];
   }
   else if (! aEdges.empty())
   {
      std::vector<int> counts(30, 0);
      for (size_t i = 0; i < aEdges.size(); ++i)
      {
         std::set<int> mats;
         int faces[64];
         size_t fCount = aAdjacency.GetEdgeFaces(aEdges[i], faces);
         for (size_t fi = 0; fi < fCount; ++fi)
         {
            int ind = faces[fi];
            ++counts[aMesh.faces[ind].materials[0]+1];
            ++counts[aMesh.faces[ind].materials[1]+1];
         }   
      }
      for (size_t i = 0; i < counts.size(); ++i)
      {
         if (counts[i] == aEdges.size() * 2)
         {
            cutRegion = i-1;
            break;
         }
      }
   }
   if (cutRegion == -999)
   {
      cout << "Cannot Detatch -- No single region shares these edges.\n";
      return;
   }

   std::map<int,int> vertexMapping;
   // duplicate faces, move offset
   for (size_t i = 0; i < aFaces.size(); ++i)
   {
      Face& oldFace = aMesh.faces[aFaces[i]];
      Face f = oldFace;
      Vec3 n = -FaceNormal(aMesh, f, cutRegion);
      for (size_t vi = 0; vi < 3; ++vi)
      {
         int v = f.vertex[vi];
         if (vertexMapping.find(v) != vertexMapping.end())
         {
            f.vertex[vi] = vertexMapping[v];
         }
         else
         {
            Vec3 vec = aMesh.verts[v];
            vec += n / 100.0;
            f.vertex[vi] = aMesh.verts.size();
            vertexMapping[v] = f.vertex[vi];
            aMesh.verts.push_back(vec);
         }
      }

      if (f.materials[0] == cutRegion)
      {
         f.materials[1] = oldFace.materials[0] = aFillMaterial;
      }
      else
      {
         f.materials[0] = oldFace.materials[1] = aFillMaterial;
      }
      aMesh.faces.push_back(f);
   }
   for (size_t i = 0; i < aEdges.size(); ++i)
   {
      int verts[] = { aEdges[i].first, aEdges[i].second };
      for (size_t vi = 0; vi < 2; ++vi)
      {
         int v = verts[vi];
         if (vertexMapping.find(v) == vertexMapping.end())
         {
            vertexMapping[v] = aMesh.verts.size();
            aMesh.verts.push_back(aMesh.verts[v]);
         }
      }
   }
   EdgeList borderEdges;
   
   std::vector<int> allFaces;
   IndexList         borderVertices;
   // Figure out which edges/vertices are on the border
   GetBorderEdges(aMesh, aFaces, borderEdges);
   GetEdgeVertices(borderEdges, borderVertices);
   
   // Find all border faces
   for (size_t i = 0; i != borderVertices.size(); ++i)
   {
      int faces[128];
      size_t fCount = aAdjacency.GetVertexFaces(borderVertices[i], faces);
      for (size_t j = 0; j < fCount; ++j)
      {
         Face::Materials& mat = aMesh.faces[faces[j]].materials;
         if (mat[0] == cutRegion || mat[1] == cutRegion)
         {
            allFaces.push_back(faces[j]);
         }
      }
   }
   make_unique(allFaces);

   // Save border faces not in selection
   std::vector<int> borderFaces;
   for (size_t i = 0; i < allFaces.size(); ++i)
   {
      int fi = allFaces[i];
      if (std::find(aFaces.begin(), aFaces.end(), fi) == aFaces.end())
      {
         borderFaces.push_back(fi);
      }
   }

   std::vector<int> deadFaces;

   // Loop through border faces, duplicate, color (cutRegion, fillRegion), Re-assign vertices using vertexMapping
   // Re-color original face (fillRegion, otherRegion), Add to delete list if new color is (fill, fill)
   for (std::vector<int>::iterator i = borderFaces.begin(); i != borderFaces.end(); ++i)
   {
      Face& oldFace = aMesh.faces[*i];
      Face  newFace = oldFace;
      for (int vi = 0; vi < 3; ++vi)
      {
         int vert = newFace.vertex[vi];
         if (vertexMapping.find(vert) != vertexMapping.end())
         {
            newFace.vertex[vi] = vertexMapping[vert];
         }
      }
      if (newFace.materials[0] == cutRegion)
      {
         newFace.materials[1] = oldFace.materials[0] = aFillMaterial;
      }
      else
      {
         newFace.materials[0] = oldFace.materials[1] = aFillMaterial;
      }
      if (oldFace.materials[0] == oldFace.materials[1])
      {
         deadFaces.push_back(*i);
      }
      aMesh.faces.push_back(newFace);
   }

   MeshAlgorithm::RemoveFaces(aMesh, FaceSelectionCriteria(deadFaces));
   //! COMPLETE!
}

void MeshAlgorithm::DetatchInterface3(Mesh&            aMesh,
                                      MeshAdjacency&   aAdjacency,
                                      CreaseStructure& aCreases,
                                      int              aPatchIndex,
                                      int              aFillMaterial)
{
   RegionCrease* crease = aCreases.mPatches[aPatchIndex];
   RegionInterface interface = crease->mInterface;
   int cutRegion = interface.first;
   int otherRegion = interface.second;

   std::map<int,int> vertexMapping;
   // duplicate faces, move offset
   for (size_t i = 0; i < crease->mFaces.size(); ++i)
   {
      Face& oldFace = aMesh.faces[crease->mFaces[i]];
      Face f = oldFace;
      Vec3 n = -FaceNormal(aMesh, f, cutRegion);
      for (size_t vi = 0; vi < 3; ++vi)
      {
         int v = f.vertex[vi];
         if (vertexMapping.find(v) != vertexMapping.end())
         {
            f.vertex[vi] = vertexMapping[v];
         }
         else
         {
            Vec3 vec = aMesh.verts[v];
            vec += n / 100.0;
            f.vertex[vi] = aMesh.verts.size();
            vertexMapping[v] = f.vertex[vi];
            aMesh.verts.push_back(vec);
         }
      }

      if (f.materials[0] == cutRegion)
      {
         f.materials[1] = oldFace.materials[0] = aFillMaterial;
      }
      else
      {
         f.materials[0] = oldFace.materials[1] = aFillMaterial;
      }
      aMesh.faces.push_back(f);
   }

   // Search loop vertices for border faces
   std::set<int> patchFaces, borderFaces;
   patchFaces.insert(crease->mFaces.begin(), crease->mFaces.end());
   
   for (size_t i = 0; i < crease->mLoops.size(); ++i)
   {
      VertexLoop& loop = *crease->mLoops[i];
      for (size_t j = 0; j < loop.mVertices.size(); ++j)
      {
         int faces[128];
         int vi = loop.mVertices[j];
         size_t faceCount = aAdjacency.GetVertexFaces(vi, faces);
         for (size_t k = 0; k < faceCount; ++k)
         {
            int fi = faces[k];
            Face& f = aMesh.faces[fi];
            // pick only faces belonging to cutRegion, but not in duplicated patch!
            if (  f.materials[0] == cutRegion && f.materials[1] != otherRegion
               || f.materials[1] == cutRegion && f.materials[0] != otherRegion)
            {
               if (patchFaces.find(fi) == patchFaces.end())
               {
                  borderFaces.insert(fi);
               }
            }
         }
      }
   }

   std::vector<int> deadFaces;

   // Loop through border faces, duplicate, color (cutRegion, fillRegion), Re-assign vertices using vertexMapping
   // Re-color original face (fillRegion, otherRegion), Add to delete list if new color is (fill, fill)
   for (std::set<int>::iterator i = borderFaces.begin(); i != borderFaces.end(); ++i)
   {
      Face& oldFace = aMesh.faces[*i];
      Face  newFace = oldFace;
      for (int vi = 0; vi < 3; ++vi)
      {
         int vert = newFace.vertex[vi];
         if (vertexMapping.find(vert) != vertexMapping.end())
         {
            newFace.vertex[vi] = vertexMapping[vert];
         }
      }
      if (newFace.materials[0] == cutRegion)
      {
         newFace.materials[1] = oldFace.materials[0] = aFillMaterial;
      }
      else
      {
         newFace.materials[0] = oldFace.materials[1] = aFillMaterial;
      }
      if (oldFace.materials[0] == oldFace.materials[1])
      {
         deadFaces.push_back(*i);
      }
      aMesh.faces.push_back(newFace);
   }
   
   MeshAlgorithm::RemoveFaces(aMesh, FaceSelectionCriteria(deadFaces));
   //! COMPLETE!
}

void MeshAlgorithm::EdgeMaterials(Mesh&             aMesh,
                                  MeshAdjacency&    aAdjacency,
                                  Edge              aEdge,
                                  std::vector<int>& aMaterials)
{
   int faces[64];
   size_t fCount = aAdjacency.GetEdgeFaces(aEdge, faces);
   assert(fCount<64);
   for (size_t i = 0; i < fCount; ++i)
   {
      Face& f = aMesh.faces[faces[i]];
      aMaterials.push_back(f.materials[0]);
      aMaterials.push_back(f.materials[1]);
   }
   make_unique(aMaterials);
}

void MeshAlgorithm::CreaseEdgePatch(Mesh&          aMesh,
                                    MeshAdjacency& aAdjacency,
                                    Edge           aEdge,
                                    EdgeList&      aEdgePatch)
{
   std::set<Edge> expanded, unexpanded;
   unexpanded.insert(aEdge);
   std::vector<int> patchMaterials, edgeMaterials;     
   EdgeMaterials(aMesh, aAdjacency, aEdge, patchMaterials);
   if (patchMaterials.size() < 3)
   {
      aEdgePatch.push_back(aEdge);
      return;
   }
   while (! unexpanded.empty())
   {
      edgeMaterials.clear();
      Edge e = *unexpanded.begin();
      unexpanded.erase(unexpanded.begin());
      expanded.insert(e);
      EdgeMaterials(aMesh, aAdjacency ,e, edgeMaterials);
      if (edgeMaterials == patchMaterials)
      {
         aEdgePatch.push_back(e);
         int v1[256], v2[256];
         int v1ct = aAdjacency.GetVertexStar(e.first, v1);
         int v2ct = aAdjacency.GetVertexStar(e.second, v2);
         assert(v1ct<32);
         assert(v2ct<32);
         std::vector<Edge> adjacent;
         for (int i = 0; i < v1ct; ++i)
         {
            adjacent.push_back(OrderEdge(v1[i], e.first));
         }
         for (int i = 0; i < v2ct; ++i)
         {
            adjacent.push_back(OrderEdge(v2[i], e.second));
         }
         for (size_t i = 0; i < adjacent.size(); ++i)
         {
            if (expanded.find(adjacent[i]) == expanded.end())
            {
               unexpanded.insert(adjacent[i]);
            }
         }
      }
   }
}
